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PROJECT SUMMARY 



The research was directed toward the development and testing 
of molecular models for ice nucleation and growth. Subsequent 
research will continue to refine the models and apply them to 
examine specific properties which influence ice processes in the 
natural atmosphere and under low gravity conditions. Specific 
short range goals were to apply and compare two distinct but com- 
plimentary techniques, the semi -empirical quantum mechanical 
approach (see Attachment I) and the central field empirical pair 
potential model (see Attachment II). The former technique was 
used to study ionic, configurational and impurity defects in the 
ice surface while the latter concentrated on studying the effects 
of kinks, ledges and vacancies on the absorption of water molecules 
from the vapor. 

In addition to static calculations of the energetics of the 
water monomer-ice surface interaction, molecular dynamics studies 
were initiated. These dynamics studies did not model the ice 
surface directly but used small aggregates of water molecules to 
examine the time evolution and temperature dependence of structure 
in hydrogen bonded systems (see Attachment I, section 2). Molecular 
dynamics simulations of collisions also provide information about 
"sticking coefficients" and the role of energy transfer in the 
absorption rrocess. 

Results of the calculations of mcnomer-ice surface interactions 
conclude that many adsorption sites exist on the ice surfaces and 
that the barriers between "bonding sites" are relatively low. Thus 
the surface of ice is exceedingly complex and these multiple 


attractive sites suggest the probability of a liquid-like layer at 
ordinary temperatures. 

Both the quantum mechanical and the empiric.' 1 potential results 
indicate that bonding on the prism face of ice 1^ is preferentially 
above lattice sites. However, for the basal surface the two tech- 
niques differ in their prediction of the most stable absorption 
sites. The quantum mechanical calculations again predict that 
bonding to a lattice site is preferred while the empirical potential 
predicts a greater stability for absorption between lattice sites. 

The source of this disagreement has not been resolved. Since the 
differences are not large, careful consideration of proton orienta- 
tions, effects of unshielued long range forces in the empirical 
model and correlation effects in the quantum model may be required. 
Then too, surface relaxation effects were not considered and these 
may alter the predictions as to the magnitude of the monomer inter- 
action with different surface sites. In addition, the molecular 
dynamics and Monte Carlo calculations are beginning to provide 
information about the effects of surface relaxation, temperature 
and collision energy transfer processes on the monomer-surface 
binding, all of which may alter the predictions of the static calcu- 
lations . 

The major results of the molecular dynamics simulations are 
i) cluster and/or surface modes play a primary role in the absorption 
process, ii) the flexibility of the hydrogen bond serves to decrease 
the energy required for structural interconversion, iii) the rapid 
distribution of added energy in a hydrogen bonded system leads to 
aggregate stability which greatly exceeds that predicted by static 
energy calculations and iv) as a result of iii) additional 


mechanisms for dissociation or evaporation processes in the nucle- 
ation kinetics should be considered. 

The results obtained during the course of this contract and 
summarized in attachments I and II show that molecular models for 
nucleation and growth of ice have been developed and have been 
successful in providing insight into these processes. Many ques- 
tions still remain and refinements of the models are needed, but 
the work to date demonstrates the viability of the molecular 
approach. 
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I. Quantum Mechanical Studies of Water Interactions with Ice 
Surfaces 

In the continuing effort to analyze the effect of microscopic 
properties of the surface on the nucleation and growth of ice, 
quantum mechanical calculations of a model ice surface interacting 
with a water monomer were undertaken. The initial calculations 
which examined several possible bonding sites on the basal and 
prism faces arc reported in ref. 1, (Plummer, 1978). The conclu- 
sions of this study were that the strength of the bonding of a 
monomer to the surface was roughly the same for both the basal and 
prism faces but that residence times and diffusion rates on the two 
faces could be substantially different. The presence of many local 
minima supgests that surface diffusion will be a primary mechanism 
for surface cluster formation and growth. 

In these calculations the surface-monomer interaction was 
modeled using a th ; rtcen molecule ice cluster with the molecules 
arranged as for a basal or prism f ace of ice lj . The technique 
used to evaluate the stability of different surfaces and configu- 
rations was a seiniempirical quantum mechanical method known by the 

i i 

acron/m CNDO/ 2 v . The ice cluster used to represent the basal 
surface is shown in Figs. 1 and 2. The oxygens represented by 

O 

filled circles lie 0.9 A below the plane. 


lig. 1. Perspective drawing of 
the basal surface of ice 1^ (BA). 
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Fig. 2. Section (^60 A 2 ) of the 
basal surface of ic 2 Iw showing 
proton configurations for an 
alternating (BA) surface. Primary 
bonding site is indicated by 1; 
dashed lines connect coplanar 
molecules and illustrate the 
hexagonal symmetry of the lattice. 



In order to better visualize the ice surface as viewed by an 
approaching water monomer, electron density plots were generated. 
Such a plot is shown in Fig. 3 for the basal face with the same 
proton configuration as shown in Figs. 1 and 2. 


Fig. 3. Electron density map of 
the basal surface of ice 1^. The 
magnitude of the density in the 
plane of the upper level o r mole- 
cules is represented by the height 
of the peaks. Total area of the 
density map is 180 A 2 ; grid points 
arc .25 A apart. Perspective is 
30° above plane and 45° clockwise 
from the orientation depicted in 
Fig. 2. 



Several possible bonding sites were investigated for the basal face. 
The results are discussed further in ref. 1 and summarized in Table I. 


TABLE 1 

Summary of Results for 

Basal (BA) Surface 

i 

-Energy Calculations 

as a 

Function of Site 

1 


Ene rgy 

| 

Frequency 

(cm' 1 ) 

Site 

( kca 1 /mole) 

Primary basal 


I 

a) 180 

-7.98 

280 I 

b) 0° 

-7.50 

! 

Second layer 

-4.05 

180 

Ring Center 

-4.45 


Bridge Site 

-6.79 

19 0 

_J 



Other neutral surfaces were also studied: a prism surface, 
two polarized basal surfaces and an optimal basal surface. The 
prism surface illustrating; the proton configuration is given as 


Fig. 4. 

Fig. 4. Perspective drawing 

p — 

b 


v 

of the 13 molecules used to 
represent the prism face of ice 

A 


/ 

1^. The primary site is indi- 
cated by the arrow. 

m 

» 

'V 



The other basal surfaces examined differ from Fig. 1 only in proton 
configuration. The polarized surfaces had all of the protons in the 
upper level pointing either up (BU) or down ( BD ) . The optimal basal 
surface had the upper protons alternating as in Fig. 1 but the 
second level protons were reoriented so as to produce the most 
stable configuration for the surface itself (BOPT). The binding 
energy of a monomer to the central site for these surfaces is given 
in Table II. 


TABLF II 


Summary of Results of Monomer -Surface Interaction Energy 
as a Function of Surface Configuration 



Site 

Binding Energy 
(kcal/mole) 

Basal 

primary site 


a) 

BA 

-7.98 

b) 

Optimal (BOPT) 

-7.46 

c) 

Proton up (BUI 

-(.62 

d) 

Proton down (BD) 

-7.13 

Pr i sm 

a) 

Primary 

-7.35 

b) 

Bridge site 

-3.88 




In addition to the neutral ice surfaces, the effect of 
including defects in the surface was examined. The first defects 
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considered were ionic defects such as ion pairs which would be 
formed by proton tunneling. The calculations attempt to simulate 
the condition of an incoming water monomer reacting with a surface 
ion site before recombination can take place. Two initial cases 
having the ions on adjacent lattice sites were considered, one 
having a positive surface site and the other a negative surface 
site. A further case modeled a second tunneling event which pro- 
duces a "separated" ion pair and a surface having both a positive 
and a negative site. The results for the ion defect surface 


studies are summarized in Table in. 


TABLE III 

Summary of Results 

for Energy Calculation for 


Monomer Interact 

ions with Defect 

Surfaces 





Rotational 


Binding Energy 

Frequency 

Barrier 

Site 

( kcal /mole ) 

(cm' 1 ) 

(kcal/molc) 

Primary basal 

- 7.98 

280 

. 50 

Defect ion pairs 




a) Positive site 

-5.07 

130 

. 20 

b) Negative site 

-24.34 

300 

.67 

c) Separated ion pair 

-37.4b 

380 

- 

Configurational Defects 




a) L-dcfect 

-7.0 

- 


b) D-^efect 

-6.8 

* 



The effect of incorporating other ions in the ice surface on 
an approaching monomer are very similar to the results for sepa- 
rated ion pairs wilh a Na impurity behaving very much like a 
hydrodium ion and a Cl impurity like a hydroxyl ion. These results 
are only indicative rather than conclusive because the ice lattice 
was not allowed to relax about the impurity ion. It might be 
expected that this relaxation might be substantial because the 
small sizes of the impurity ions. 
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The other types of defects modeled were the so-called config- 
urational defects or D- and L- defects. These types of defects 
occur when a molecular rotation becomes "frozen in" by the incor- 
porating of additional lattice members. In such defects, there are 
either two (D-type) or no (L-type) protons lying along the normal 
hydrogen bond. Since such defects occur between oxygen layers they 
are not as close to an incoming monomer and therefore the initial 
effect on the interaction is considerably less than for either a 
polar surface or ion defect. However the influence of the config- 
urational defect would be expected to be greater for secondary 
bending sites and for surface diffusion because of the lattice 
relaxation about the defect. These results are also given in 
Table III. 
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I.. Molecular Dynamics Studies 

The quantum mechanical calculations described in the previous 
paragraphs served to provide much information about the monomer - 
ice surface but at an effective temperature of zero kelvin. To 
examine the effects of non-zero temperatures and also to examine 
t ime -dependent phenomena, a series of molecular dynamics studies 
were undertaken. Because of the complexity of the ice problem, 
these studies were carried out for small clusters of water mole- 
cules. Using clusters rather than a model for bulk ice, the 
importance of different energy transfer processes, molecular 
motions, and molecular orientations, including proton positions 
could be examined. Since almost all cluster molecules weie also 
"surface molecules" the ef r ects being studied are more related to 
surface than bulk processes . 

In the first series of studies (see the attached preprint"’ 

and the dissertation of T. S. Chen) the structure anJ stability of 

a twenty molecule clatnrate was investigated. The interaction 

4 

potential used was a centr°l field pair potential . The equations 
are given below. 

l'Ui.SJ8 26758.2 9 , 9 

Voq(r) ♦ - -25 exp [ - 4(r-3.4)* ) - .25 exp[ -1 .5(R-n .5)^ ] 


y (R) . . 72.269, 6.23003 10 . 0 

r R 9. 19912 40(r 1.05)] 1 ♦ ; xpI 5.4s505(r-2 .2)) 


V mm (r ) - 36 .1 19 5 ♦ IS 17 t xr ( -7. 621 77(r-1. 45251)^1 

R 1 ♦ lxp[ 40(r-2.05) J 
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The molecules were initially arranged on the vertices of a regular 
dodecahedron (Fig. 5). 

Fig. 5. Symmetric proton 
configuration for the 
clathrate (pentagonal 
dodecahedron) . 


In this cluster all of the water molecules participate in three 
bonds, as in a fully bonded surface layer. The stability of the 
structure was examined as a function of temperature, proton 
arrangement and method of energy addition to the system. It was 
found that the cluster was remarkably stable to dissociation even 
under conditions of extreme heating. The structure or hydrogen- 
bonding pattern was disrupted above 250K but an evaporation event 
did not occur until much higher temperatures. The excess energy 
was absorbed by disruption of the hydrogen bond structure (melt- 
ing) with bonds continually breaking and reforming. When dissoci- 
ation finally occurred, it was not through the evaporation of a 
single molecule but via the loss of clusters containing two or more 
molecules. Upon examination of the mechanisms for energy trans- 
mittal through the cluster, the molecular dynamics shows that the 
energy does not build up in the low frequency modes, which would 
ead to dissociatior but rather is transmitted by the librational 
modes throughout the cluster. This transfer of energy takes place 
extremely rapidly (< .02 psec) and appears to account for the 
remarkable stability of the hydrogen bonded system. The method of 
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break-up of the cluster suggests that usual nucleation kinetics 
model which allows for decay by only one molecule should be re- 
examined if this trend is observed for other cluster sizes as well. 

To begin to model the nucleation and growth processes a 
second set of studies of collision and scattering simulation between 
the clathrate and a monomer were undertaken. Four types of colli- 
sion processes were modeled. The random proton structure was used 
and the initial cluster temperature was ~ 225K. 

1) A virtual collision was simulated by giving one cluster mole- 
cule extra translational energy--of the order of 2.5 kcal/mole. 

In this process the injected energy was directed toward the center 
of mass cf the cluster, and the energy redistribution mechanisms 
were studied. The sudden transfer of energy to a single molecule 
did not, as has been speculated, disrupt the cluster structure sig- 
nificantly nor did it result in the loss of a molecule. The energy 
was rapidly transferred to adjacent molecules and was absorbed by 
the cluster. Subsequent observations revealed a radial distribution 
function relatively unchanged from its initial form before the 
virtual collision. There was no evidence for any "build-up" of 
energy which could result in the evaporation of a molecule. 

2) A virtual evaporation was simulated by instantaneously removing 
one molecule from the cluster. The remaining molecules were left 
unaltered. Surprisingly, the cage-like structure did not collapse 
and there was no extensive rearrangement of the remaining bonds. 

The cage expanded somewhat, however, and the radial distribution 
showed loss of definition in the third and subsequent peaks. Much 
of the nearest neignbor and f ive-membered ring structure of the 
original cluster remained intact. 
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3) An additional collision was simulated by introducing a water 
molecule as a projectile with an average kinetic energy for the 
temperature under study, and with impact parameters of 0, 5 and 
10 A. (The cluster radius was about 5 A.) In the first of these 
collisions the incoming monomer was allowed to transfer energy 
internally but not to the cluster, which was held rigid. Colli- 
sion with a rigid cluster attempted to mimic a collision with a 
very cold cluster and also to examine the ability of the projectile 
to accommodate to the cluster in an attempt to maximize the inter- 
action. The results of these collisions arc shown in Figs, b, 7 
and 8. As the attraction between the cluster and the monomer 
increased there was a corresponding increase in the monomer libra- 
tional energy and decrease in monomer translational energy--ref loo- 
ting accommodation of the monomer to the cluster. In all cases, 
the monomer eventually escaped. At an impact parameter of 0 A 
(Fig. b) the monomer bounced once before escaping. 

Fig. b. Collision sequence of 
a monomer with a rigid 20 molecule 
cluster. The number alongside the 
trajectory ir. elapsed time in 
units of 1,000 At. The impact 
parameter, b, is OA. 



At 5 A (Fig. 7) the monomer scatters off t he cluster. 



Fig. 7. Collision of a 
monomer with a rigid 
c lathrate (b = 5 A) . 
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At 10 A an orbiting collision was observed (Fig. S'). 

Fig. 8. Collision of a water 
monomer with a rigid cluster 
(b = 10 A). The lines along 
the trajectory indicate the 
changing monomer orientations. 


In the latter case the interaction energy was large-- in excess of 
that required for bond formation. However, large oscillations in 
the librational energy introduced reorientation to .» weaker inter- 
action configuration and the monomer subsequent lv escaped. 

4") This study was the same as case 3 but now included the full 
dynamics of the target cluster as well as the projectile. Again 
initial impact parameters of 0 A, 5 A or 111 A were used. lor all 
impact parameters, the impinging monomer was incorporated into the 
cluster-- indicat ing that energy transfer from the monomer to the 
cluster molecules was an efficient process. The interaction 
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energies (between the monomer and cluster) ranged from 7 to 10 
kcal/mole, or 1.5 to 3 hydrogen bonds. The cluster structure was 
altered under the influence of the monomer even when the interaction 
energy was small (less than one kcal/mole) . In the addition process, 
the cluster absorbs about 10 to 15% of the initial monomer energy. 


Figures 9 and 10 indicate two "snapshot” views of the monomer- 
cluster absorption process for the 0 A impact parameter case. 


lw - 215 X 

Impact - OA 
1 1 me - 1254 t 
l - 0.0 


* 



Fig. 9. "Snapshot" drawing of 
the monomer and cluster configu- 
rations. The impact parameter 
is 0 A, t=0. 



Fig. 10. Subsequent view of 
the monomer and cluster con- 
figurations for the case with 
impact parameter equal to 0 A, 
t=l psec. 


For all impact parameters the results indicated that once a bond was 
formed (when the average interaction potential exceeds 4 kcal/mole) 
the monomer reorients within about 100 intramolecular ibrations to 
form a second bond with the cluster. 

A third series of molecular dynamic studies was made on even 
smaller water aggregates. Clusters containing five water molecules 
were used to examine the rate and temperature dependence of struc- 
ture interconversion. First a series of different hydrogen bonding 
clusters were cooled to determine which geometries were the most 
stable at low temperatures (1-5 K) . These structures are shown in 




Fig. 11. 


Fig. 11. Optimization of 
various water pentamer 
structures. On the left 
side of the arrow are 
initial configurations. 

On the right side of the 
arrow are optimized 
configurations. 
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Then a series of heating and equilibrium runs were undertaken as 


indicated in the flow chart in Fig. 12. 


Fig. 12. Flow chart for 
structure stabili ty 
simulations. 
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It was found that the cage structure was the most stable at low 
temperatures. However, the barriers between the "stable structures" 
were sufficiently high that no interconversions occurred below about 
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70 K. Thus we conclude that at low temperatures, the configuration 
space is not ergodic and all such "stable" structures must be 
included in statistical studies. 

As expected, at the low temperatures, the energy of the hydro- 
gen bonds dominate the free energy for the aggregate. As the 
temperature is increased, the entropy term becomes increasingly 
more important resulting in the transition from the cage to the 
fused ring to the ring structure. The conversion process is indi- 
cated in Fig. 13 by a plot of the radial distribution function as 



Fig. 13. The plots show the change in the 0-0 distribution as the 
case structure is heated and undergoes conversion first to a fused 
ring and then to a ring. 


Fig. 14 shows the energy vs. temperature for two structures, C-, and 
ring. 
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Fig. 14. The total energy of 
the ring pentamer as a function 
of temperature is plotted. The 
structure or bonding pattern 
remains essentially unchanged 
from 2 K to 280 K. Above 290 K 
structures having fewer total 
hydrogen bonds or "dangling" 
molecules are observed with 
increasing frequency. The 
change in slope around 300 K 
corresponds to fragmentation of 
the ring and possibly can be 
considered as a transition from 
dense to a gas phase. The 
other points on the curve show 
the transition of another 
pentamer structure, C2, to the 
ring at approximately 100 K. 

The C2 structure formed on 
cooling a tetrahedral pentamer. 


Above 29UK structures having fewer total bvdrogen bonds or "dangling" 
molecules are observed with increasing frequency. The change in 
slope around 300K corresponds to fragmentation of the ring and 
possibly can be considered as a transition from a dense to a gas 
phase. For the ring the fraction of hydrogen bonds having a specific 
length at 20K and 297K are shown in Fig. 15. From this and the 
frequency spectra at the same temperatures (Fig. 10, 17) it is seen 
that although the average structure is still a ring at the higher 
temperatures, the hydrogen bonds are becoming very "fluid" and an 
increasing amount of energy is contained in the librational motions. 
As was seen before the case of the larger clusters, a buildup of 
energy in the librational modes leads to a disruption of the hydro- 
gen bonding structure. 
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Fig. 15. Plot of the 
fraction of 0-0 separation 
at a given distance for 
the riig pentamer at 20 K 
and 297 K. 



Fig. 16. Frequencies 
associated with group 
cluster motions of the 
pentamer ring at 20 K. 


Fig. 17. Same as 16 but 
at 297. Note the change 
of scale on the vertical 
axis. 
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Detailed analysis of these motions at higher temperatures has shown 
how readily certain structural changes and proton switching reac- 
tions can occur, even when the energy present is less than that 
required to break a hydrogen bond. This is illustrated in Fig. 18 
which shows the librational motion which results in exchanging the 
hydrogen participating in the hydrogen bond with the non bonded 
hydrogen and in Fig. 19 which shows the energy components for the 


exchange . 


Fig. 18. The sequence of 
orientations observed for 
proton switching reaction 
in a pentamer ring at 
- 220 K. 



Fig. 19. Plot of energy vs. 
time for the electrostatic 
and total potential energy 
during a proton switching 
reaction . 
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Thus the energy required for the switch was 3 kcal/mole rather 
than 7 kcal/mole. Such switching reactions and cooperative effects 
serve to substantially alter the energy requirements for structure 
interconversions and hence for diffusional motion and surface 
aggregation . 

In a further attempt to assess the ease of structural changes 
during the growth of an aggregate, collision studies between a 
monomer and/or pentamer were made. Typical results are shown in 
Figs. 20 and 21. 


Fig. 20. Collision sequence 
of monomer with pentamer. 
Impact parameter 6 A. 


Fig. 21. Collision sequence 
of monomer with ring 
pentamer. Initial monomer 
position 6 A above ring. 
Impact parameter 0 A. 
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The incoming monomer initially had an energy equal to the average 
energy of a cluster molecule. Trajectories examined were i) in 
the plane of the cluster at impact parameter OA and ii) b A and 
iii) perpendicular to the plane of the cluster with impact param- 
eter of OA. Here again the monomer was captured and quickly formed 
a second hydrogen bond with the cluster. Ring enlargement reactions 
were observed with the perpendicular trajectory indicating the 
transition from a 5 to a 6 member ring readily occurs. 
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III. Summary of Results 

These studies have provided the techniques for studying the 
molecular dynamics of a water molecule colliding with an ice sur- 
face . Such studies can examine the kinetic factors needed to 
evaluate the nucleation rate expressions. They can also be used to 
examine the structure of the surface as a function of time. It is 
V;r>3d that they can quantify the idea of a "liquid-like" layer on 
the ice surface: at what temperature does it begin to form, how 

disordered is the surface, through how many layers does the dis- 
order (liquid layer) persist? In addition, the molecular dynamics 
has the potential of evaluating the importance of the binding sites 
identified by the static calculations to the formation of critical 
embryos and in the crystal growth process. 

These studies of the specific features of the interactions of 
water molecules with ice surfaces and clusters have provided the 
following : 

1) There are many sites on the basal and prism faces of ice at 
which the monomer is at least quasi-bound. 

2) The magnitude of the bonding on the basal and prism sites is the 
same but the barriers for moving between the sites appear to be lower 
for the basal face. This would suggest a difference in activity for 
th' two faces which could well he temperature dependent. 

3) The molecular dyi v.-ics of vapor phase clusters has shown that a 
variety of configurations for a given size cluster can co-exist. 

It is expected that this effect would be enhanced for embryos on a 
surface . 

4) The dynamic studies of monomer-cluster scattering suggest the 
possibility of larger absorption coefficients, of the order of one. 
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5) The monomer-cluster scattering experiments have also shown that 
cluster structure can be substantially changed (for example from a 
five to a six member ring - ) and that only a small percent of the 
monomer energy needs to be transferred to the cluster for absorption 
to take place. 

Thus this work with the moleculir model has shown both the 
feasibility and importance of the approach, has provided insights 
as to the importance of specific effects and has provided potentials 
and algorithms for use in the further studies of the formation and 
growth of ice. 
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ference on Chemical Dynamics, Univ. of Texas- Aust in , March, 1981. 
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Session Chairman: Symposium on Nucleation Theory and Experiment, 

ACS, March, 1977. 

Session Rapportuer: International Conference on Physics and 

Chemistry of Ice, England, September, 1977. 

Co-Organizer: Nucleation Symposium, Spring Meeting of American 

Chemical Society, 1977. 

CAUSE Proposal Review Panel, Washington, D.C., January, 1978. 

Invited participant: Symposium in honor of Robert Mulliken, Chicago, 

IL. October, 1978. 

Session Chairman: Zettlemoyer Symposium on Nucleation, 53rd Col- 

loid and Surface Science Symposium, June, 1979. 

NSF Women in Science, Proposal Review Panel, Washington, D. C., 
February, 1979, panel leader. 

Co-Organizer: Zettlemoyer Symposium on Nucleation; Symposium on 

Light Scattering; Symposium on Aerosols, Annual meeting of Col- 
loid and Interface Science Division of ACS, June, 1979. 

Invited participant: Workshop on Stochastic Molecular Dynamics, 

sponsored by NRCC , Woods Hole, Mass., July, 1979. 

Invited participant: Workshop on Computational Methods for 

Molecular Structure Determination: Theory and Technique, sponsored 

by NRCC and QCPE, Indiana Univ., August, 1979. 

Conference on Atmospheric Water Vapor; Vail, CO, September, 1979. 

Invited participant: Horizons in Hydrogen Bond Research, Sanga- 

Saby, Sweden, August, 1980. 

NSF Women in Science, Proposal Review Panel, Washington, L'. C., 
February, 1981, panel leader. 

Conference Organizer, VI International Symposium on Physics and 
Chemistry of Ice, to be held 2-6 August, 1982, Rolla, MO. 

Editorial Advisory Board: Journal of Colloid and Interface Science, 

1980-1985. 
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I. STUDIES OF THE WATER MOLECULE BONDING TO ICE SURFACES 

Our attention was focused on the microscopic aspects of ice crystal growth 
habit and ice nucleation because of the essential role these processes play in 
cloud microphysics and because of the absence of reliable macroscopic models. 

At a 1978 Workshop on Ice and Snow Crystal Growth, Weickmann^ reviewed the great 
variety of ice crystal growth habits and commented on the lack of a unified 
explanation which could incorporate: the role of a transition (or liquid like) 

layer on the ice surfaces; the role of supersaturation with respect to ice at 
water saturation values as a function of temperature; the role of two dimensional 
nucleation processes on the basal and prism faces of ice. 

During the grant period we have developed a program to examine the basic 
interaction of the water monomer with a variety of ice surfaces. These studies 
were intended to give a molecular description of the bonding and diffusion of 
H^O on ice substrates with experimental values for lattice constants and struc- 
ture. A companion study is being made on the interaction of the water molecule 
(and clusters of water molecules) on the basal and prism face of silver iodide 
The latter work is funded by grant from the National Science Foundation. 

In the first set of studies, molecular arrays of point charge water mole- 
cules are generated to represent the basal and prism faces of hexagonal ice. 

The adsorbed water monomer and the water molecules in the ice surfaces interact 
via the H ? 0-H20 potentials of Stillinger and Rahman"*. 
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where H and 0 refer to the hydrogen and oxygen components of the water molecule, 
and r is the separation distance between the atoms in the subscripts of V. These 


1 



2 


potentials have been tested for application to liquid water systems. Further 
discussion of the potentials can be found in the literature^. 

The prism (basal) faces of ice are represented by 100 (50) water molecules 
arranged in four (two) layers of 25 molecules each in the hexagonal ice 1^ 
structure. Lattice parameters are a ■ A. 519 A and c ■ 7.330 A, and each water 
molecule (with oxygen on a lattice point) is surrounded tetrahedrally by four 
water molecules 2.76 A distant. All H^O are rigid with H-O-H bond angle equal 
to 109.5° and OH distance euql to 0.9584 A. The hydrogens in the model sub- 
strate are arranged 90 that the top surface layer has water molecules with dipole 
moment projection alternately into and out of the surface plane. Basal ice faces 
with positive and negative surface polarization have also been modelled^. The 
ice substrate models are finite in size to avoid complicated Ewald summation** 
and to prevent divergences associated with periodic polarized unit cells. The 

adsorbed water molecule interacts with all the ice substracte molecules, but is 

2 

constrained with its center of mass projection over a small region about 30 A 
in area in the center of the array. Since we are interested in the qualitative 
features of the interactions and not exact magnitudes of binding energies the 
finite ice model serves adequately. Larger arrays have been examined and found 
to produce the same qualitative results. 

A four layer ledge on the ice basal face is represented by 120 water mole- 
cules arranged in four layers of 15 molecules each and based on two layers of 25 
molecules described above. The upper and lower terrace (basal faces) are nega- 
tively polarized and the prism face oolarization is the same as that described 
above. Proton configurations for all of the ice substrate models can be seen in 
the figures displaying the results. 

The maximal binding energy surfaces are generated by minimizing the potential 
energy of the adsorbed water molecule at a set of grid points about 0.25 A apart 
and covering the "adsorption" region in the center of the arrays. The center of 
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mass of the water monomer is held with fixed (x,y) projection onto the grid point, 
and the center of mass height above the substrate, z, and the monomer Euler angles 
(a, 8,^) are varied until the interaction potential energy is minimized. The mini- 

g 

mizatlon is done with subroutine VA04A of the Harwell subroutine library. 

Results for the water monomer on the prism (basal) face are shown in Fig. 2a 
(b) in the form of contour density plots. Shading In the plots indicates regions 
of higher binding energy for the adsorbed water monomer. Circles on the plots 
indicate location and proton arrangement of H^O molecules in the first and second 
layer of the ice substrate. The largest maximal binding energy site on the two 
faces is on the prism face (10.2 kcal/mole) and is located in the center of the 
six membered ring shown in the upper part of the figure. The smallest binding 
energy is also on the prism face (2.7 kcal/mole) — directly over molecules in the 
second layer. An interesting view of all the four ice substrates studied is 
shown in Fig. 1 where the minimal potential energy of the adsorbed water monomer 
is plotted three dimensionally. Wells in the potential surface indicate regions 
of large binding energy. Small dots indicate positions of ice substrate molecules 
in the first layer. In order to examine the effect of extended regions of larger 
or smaller binding energies, an average maximal binding energy is calculated for 
each model substrate by weighting the value at each grid point with a Boltzman 
factor. The averages for the prism and basal faces are 8.9 and 8.3 kcal/mole, 
respectively. For the negatively and positively polarized basal faces the results 
are 7.0 and 6.7 kcal/mole, ’•espec tively . Thus the prism face of ice appears to be 
more attractive to the adsorbed water molecule than the basal faces. The energy 
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Fig. 1. Perspective views of the minimal potential energy surfaces of the adsorbed 
water monomer on the basal and prism faces of ice. Sample H 2 O configurations with 
maximal binding energies (in kc?l/mole) are indicated. Small dots on the contours 
show positions on first layer ice substrate water molecules. 







Fig. 2. The maximal binding energy of the water molecules above sites on 
the (a) prism face and (b) unpolarized basal face of ice Ij,. The H 2 O 
molecules in the z ■ 0 plane are indicated by solid circles; z - -0.9 A 
plane molecules are Indicated by dashed circles. The smallest binding 
energy contours are 5.0 kcal/mole lor the prism face and 6.5 kcal/mole 
lor the basal face. Neighboring contours differ by 0.5 kcal/mole. The 
largest maximal binding energies ire 9.5 and 10.2 kcal/mole for the basal 
and prism faces, respectively. 

Poor 


Fig. 3. The maximal binding energy 
contours lor the water monomer above 
sites on the basal face ledge. Solid 
circles on and below the dashed line 
of the ledge are molecules in the 
first laver oi the upper terrace. The 
lower terrace is at z = -7.33 A, and 
located in the upper portion of the 
1 I gure. 
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barriers to diffusion from grid point to grid point on the maximal binding energy 
surfaces are 3.0 and 2.5 kcal/mole for the prism and basal surfaces, respectively. 
Diffusion barriers for the basal faces appear to be Independent of the surface 
polarization. 

Maximal binding energy for the water monomer adsorbed on the ice basal face 
ledge is shown in Fig. 3. The largest maximal binding energy site is at a posi- 
tion about halfway up the step — in a configuration close to that observed for 
optimal binding on the prism face without the ledge. The maximal binding energ' 
is 10.2 kcal/mole and only slightly larger than that found for the optimal binding 

configuration on the smooth prism face. The potential energy surface and some 

9 

optimal bonding configurations are shown In Figs. 4a and 4b. 

The main features of the maximal binding energy contours are: (1) potential 

energy "troughs" which appear to pass around and between first layer ice sub- 
strate molecules; (2) a preference for water monomer adsorption on the prism face 
and above basal face ledges; (3) smaller diffusion barriers on the basal faces; 

and (4) a complexity oi preferred adsorption sites with the adsorbed greatly 

affected by neighboring surface proton arrangement. The results of this model 
present a somewhat different view of H^O attaching to an ice substrate, and point 
out the difficulties in modelling the H^O - ice substrate system with the assump- 
tion that molecules adsorb directly above substrate molecules in bulk ice config- 
urations. This calculation indicates that a more '’ompl i< •• ted situation could 
prevail — particularly at low monomer concentrations — and hints at the very real 
possibility of a quasi-liquid like layer on the surface of ice. It is also noted 
that two dimensional ice clusters forming on the rigid (low temperature) ice 

surface could assume considerable distortion along the cluster periphery. Since 

it is the edge energy that determines the energy of formation of the cluster on 
the substrate, two dimensional nucleation (with some modifications) could once 
again be a possibility. This was suggested some Lime ago by Bartley^. Future 
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work will Investigate binding energies of adsorbed ice clusters and attempt to 
determine critical cluster sizes . 11 
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II. MOLECULAR DYNAMICS STUDIES OF THE WATER MONOMER AND SMALL WATER CLUSTERS 

ON THE BASAL FACE OF ICE 

These initial studies of the water monomer and a six molecule water cluster 
on the basal face of ice were begun in order to examine the diffusion rates and 
the stability of small water aggregates on the ice surface. The basal face ice 
surface (neutrally polarized) is constructed as in reference 9 , and the water 
molecules interact via the Stillinger potentials^ * ^Molecular dynamics is used to 
follow the development of the water monomer and the six molecule cluster. In 
Fig. 5 is shown a six-molecule water cluster in the process of breaking up 

after having been placed on the ice surface in an extended bulk ice lattice con- 
figuration. The six membered ring cluster was placed on the surface with zero 
kinetic energy, but at a position slightly offset from a relative minimum in the 
interaction potential energy. The displacement from the relative minimum was 
small — amounting to about 3% of the total binding energy of the cluster to the 
surface. This deviation in potential is about 1 kT at a temperature of 260 K, 
and is of the order of thermal fluctuations of the cluster. Figure 5 shows the 
water molecules diffusing to sites on the ice surface which are optimal adsorption 

sites. These optimal adsorption sites were studied in the static studies of the 

3 9 

water monomer interacting with the ice surfaces.’ We are presently extending 
these studies to allow the ice surface molecules to move and to study the exchange 
of energy between the cluster molecules and the ice surface. The results to date 
indicate that the six membered ring (ice-like) structure is not stable on the ice 
surface for times small compared to the time between collisions with adsorbed 
surface monomers. A study of the optimal binding energy surfaces for the adsorbed 
water monomer has indicated that molecules on the periphery of the surface clus- 
ters will find stronger binding to the surface at "interstitial" sites in dis- 
torted bonding configurations.* Since the six-membere-’ ring cluster is composed 
entirely of "edge" molecules, it is not surprising that the structure is unstable 
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Fig. 5. A 9 ix-membered ring water cluster adsorbed on the basal face of ice. 

O circles represent water molecules in the ice substrate; @ molecules in 
the cluster at time t*0; @ molecules in the "cluster" at t 3 10~ 13 sec. 
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on the Ice .urface. We are presently developing computer programs to look at 
larger ice-like clusters on the ice substrates to investigate the dependence of 


stability on cluster size. 
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III. STUDIES OF A PERIODIC ICE STRUCTURE USING THE STILLINGER POTENTIALS 

These studies were undertaken in order to examine the equilibrium ice struc- 
ture resulting from the Stillinger potentials^ For this investigation, a 192 
molecule ice sample was generated which has zero dipole, quadrupole and octupole 
moments. 13 The ice sample also has a hydrogen arrangement which allows the 192 
molecules to form a "unit cell" for application of periodic boundary conditions. 

The latter conditions are necessary for simulations of bulk ice. 

The ice sample was allowed to relax to its equilibrium configuration (at 
0°K) and the corresponding lattice constants for the hexagonal ice structure were 
recorded. The internal 0-H distances and H-O-H bond angles were also allowed to 
vary. These parameters were then used to generate a new unit cell for the 
periodic ice structure. 13 In the next phase of this study the unit cell in the 
minimal energy ice structure was allowed to equilibriate at 20°K using a 
Metropolis Monte Carlo method l2 . Figure 6 shows one layer of the unit cell after 
50 thousand steps. The darkened molecules represent the Initial configuration. 
There is considerable relaxation at the surface of this large ice sample — 
reflecting the deficiency of nearest neighbors in the surface molecules. 

Present studies are examining the relaxation of the periodic ice structure 
using he Monte Carlo method at 20°K. The results indicate that the initial 
minimization, which produced values 2.78 A, 0.97 A, and 101° for the 0-0 separation 
distance in the hexagonal ice, the 0-H distance internal to the H 2 O molecule, and 
the H-O-H angle internal to the water molecule, respectively, is consistent with 
a full Monte Carlo equilibration at 20°K. That is, the Stillinger potentials 
support a hexagonal ice structure at this temperature with the above constants 
generating the lattice constants and the internal H 2 O parameters. Subsequent 
work will investigate the temperature range in which this structure is stable 
and examine the relaxed surface of the sample cleaved along basal and prism 


planes . 




Fig. 6. One layer of the basic cell of the ice sample after 50 thousand Monte 
Carlo steps at 20°K. The darkened molecules indicate the initial configurations. 
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IV. COMPANION STUDIES OF WATER ON AN ICE NUCLEATING AGENT 

Under National Science Foundation Grant3 ATM77-12614 and ATM80-15790 a com- 
panion study of ice nucleation on silver iodide has been carried out. The 
objectives have been to determine critical cluster sizes, and nucleation rates 
for water and ice on hexagonal Agl, to examine the properties of water layers on 
the same substrate and to determine which substrate features promote ice forma- 
tion. For results on these studies we refer the reader to the attached copies 
of references 2, 3 and 4. The studies primarily indicate a preference for mono- 
layer water clusters containing five and six membered rings. Diffusion rates on 
the basal Agl faces were found to be larger than on the prism face and binding 
to the substrate was preferred on the prism face — at sites not directly above 
first or second layer substrate atoms. 

An example of a two layer water cluster on the Agl basal face is shown in 
Fig. 7. The upper part of the figure shows the layering (viewed paral lei to the basal 
face) and the lower portion shows the predominance of five and six membered 
rings centered on the exposed iodine atoms. This figure was generated from a 
Metropolis Monte Carlo simulation of 88 water molecules adsorbed on Agl at 265°K. 
Similar studies of water layers on the ice substrates are being developed — using 
the periodic ice sample described in the previous section. Critical cluster sizes 
on the Agl basal and prism faces are also being studied — using the technique 


described in reference 11. 
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Fig. 7. Snapshot of a two layer water system on the basal face of silver 
iodide. Top view is parallel to the substrate; bottom view is looking down 
on the basal face. Evidence of six-membered ring structure is .;oen in the 
top layer of water molecules. 
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"The Water Monomer on the Basal Plane of Ice Ij,: an Effective Pair, Central 
Force Potential Model of Static Interaction", J. Kiefer and B. Hale, J. Chem. 
Phys. 67, 3206 (1977). 

"Studies of Water on Agl Surfaces: an Effective Pair Potential Model", J. Kiefer 
and B. Hale, J. Chem. Phys. 73, 923 (1980). 

"Theoretical Studies of Water Adsorbed on Silver Iodide" B. Hale, J. Kiefer, S. 
Terrazas, and R. Ward, J. Phys. Chem. 89, 1473 (1980). 

"The Water Monomer on Surfaces of Ice 1^: Comparison of the Basal and Prism 
Faces", B. Hale, J. Kiefer, and C. Ward, J. Chem. Phys., _75, (August 15, 1981). 

"A Monte Carlo Method for Approximating Critical Cluster Size in the Nucleation 
of Model Systems", B. Hale and R. Ward, submitted to J. Statistical Physics. 
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R. Ward, J. Holdman and B. Hale, submitted to J. Chem. Phys. 
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"A Lattice Gas Model for Adsorption on Agl", S. Terrazas and B. Halt*, Missouri 
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WMO/IAMAP Scientific Conference on Weather Modification, Tashkent, U.S.S.R., 
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! 

I Ice Nucleation Workshop, Laramie, Wyoming, June, 1975. 

International Cloud Physics Conference, Boulder, Colorado, July, 1976. 

Second '/MO Scientific Conference on Weather Modification, Boulder, Colorado, 
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I 

Nucleation Symposium, American Chem. Soc . Meeting, New Orleans, March, 197.'. 

p 

Midwest Soli'.. State Conference, Purdue University, October, 1977. 

NSF Review Panel for Proposals on Women in Science, January, 1977. 

Ice Crystal Growth Habit Workshop, Boulder, Colorado, February, 1978. 

American Chemical Society Meeting, Anaheim, California, March, 1978. 

Gordon Conference on "Chemistry at Interfaces", Meriden, New Hampshire, July, 
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Visit to Institut for Experimental Physics, University of Vienna, September, 
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Gordon Conference on Water and Aqueous Solutions, Plymouth, NH, July 1980. 
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X. SUMMARY OF PROGRESS 

During the grant period we have developed molecular models for the water 
molecule's environment on the positively, negatively and neutrally polarized 

Q 

ice I basal faces/ on the neutrally polarized ice I prism face, and or. a 

negatively polarized ice basal face with a four-molecular layer ledge. Q The 

water monomer was found to adsorb more strongly on the prism face than on the 

basal faces of ice — at sites not directly above molecules in the first or 

second layer of the (rigid) ice I. substrate. The basal face four layer ledge 

n 

served to expose favorable sites on the prism face and to promote attachment of 
the H^O at the base of the step. Diffusion rates fr- the monomer wen found 
to bt slightly larger on the basal faces 7 » 'h’hese results, together with the 
multiplicity of relative maxima found in the binding energy for H^O adsorbed on 
ice, indicate a rather complex environment for the adsorbed water. The latter 
is consistent with the existence of a liquid-like layer on the ice surfaces. 1 
Molecular dynamics studies of the decay of a six membered ring ice structure on 
the (rigid) ice basal further strengthen this hypothesis. A 192 molecule 

periodic ice structure was generated which has zero dipole, quadrupole and 

1 ? 

octupole moments. The sample was allowed to relax to an equilibrium hexagonal 
ice lattice using the Stl lliager central iorce potentials. The optimal lattice 
constant values were determined, together with optimal values for the H^O 
internal 0-H distance and H-O-H bond angle. 11 This energy minimized periodic ice 
structure is presently being studied via Metropolis Monte Carlo methods to deter- 
mine stability as a function of temperature. The surface of this ice structure 
will be applied to studies of the stability and structure of basal and prism icesurlace 
layers and subsequently to studies of ice crystal growth habit and nucleation 
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